In all of the foregoing, we have assumed a fixed value of the dispersion parameter \(k\), of the negative binomial measurement model. Now relax this assumption: estimate \(k\) along with the other parameters. How much is the fit improved? How do the MLE values of the other parameters change?
First, we set up the problem, loading packages, initializing a parallel environment, constructing our SIR model, and identifying the putative MLE. The pomp constructed by the source command below is named measSIR. We take for our putative MLE the point in our database with the highest measured likelihood.
library(pomp)
library(tidyverse)
library(doParallel)
library(doRNG)
registerDoParallel()
source("https://kingaa.github.io/sbied/pfilter/model.R")
read_csv("measles_params.csv") %>%
filter(!is.na(loglik), loglik.se<1) %>%
filter(loglik==max(loglik)) %>%
select(-loglik,-loglik.se) -> coef(measSIR)
The putative MLE is
| Beta | 3.79e+00 |
| mu_IR | 2.00e+00 |
| rho | 5.79e-02 |
| k | 1.00e+01 |
| eta | 5.88e-01 |
| N | 3.80e+04 |
We’ll begin, as in the Lesson, by performing a local search with only \(N\) fixed. The purpose of this is to evaluate the performance of the search algorithm, so that we can tune the algorithm’s parameters.
Windows users, beware! The following triggers a C snippet compilation within a parallel block, which will not succeed on many Windows machines. See this discussion for an explanation of how to circumvent this problem.
registerDoRNG(482942941)
foreach(i=1:8,.combine=c) %dopar% {
library(pomp)
library(tidyverse)
measSIR %>%
mif2(
Np=2000, Nmif=40,
cooling.fraction.50=0.5,
rw.sd=rw.sd(Beta=0.02, rho=0.02, mu_IR=0.02, k=0.02, eta=ivp(0.05)),
partrans=parameter_trans(log=c("Beta","k","mu_IR"),logit=c("rho","eta")),
paramnames=c("Beta","rho","k","eta","mu_IR")
)
} -> mifs_local
We examine the traces of the IF2 runs:
mifs_local %>%
traces() %>%
melt() %>%
filter(variable != "N") %>%
ggplot(aes(x=iteration,y=value,group=L1,color=factor(L1)))+
geom_line()+
guides(color="none")+
facet_wrap(~variable,scales="free_y")
IF2 seems to be doing a reasonable job of increasing the likelihood, though there is some way yet to go. We note that \(k\) is immediately decreased, giving the model more measurement error, with which to explain discrepancies between the process model and the data. Of course, to maximize the likelihood, \(k\) should be as large as possible. It will be interesting to see what happens to the other parameters as IF2 tries to increase \(k\).
We now evaluate the log likelihood at each of the resulting points.
registerDoRNG(908222057)
foreach(mf=mifs_local,.combine=rbind) %dopar% {
library(pomp)
library(tidyverse)
evals <- replicate(10, logLik(pfilter(mf,Np=2000)))
ll <- logmeanexp(evals,se=TRUE)
mf %>% coef() %>% bind_rows() %>%
bind_cols(loglik=ll[1],loglik.se=ll[2])
} -> results
The above calculations took a total of 12Â s on 20Â cpus. What do they show us?
pairs(~loglik+k+Beta+eta+rho+mu_IR,data=results,pch=16)
As usual, we store these results in a database.
read_csv("measles_params.csv") %>%
bind_rows(results) %>%
arrange(-loglik) %>%
write_csv("fitall_params.csv")
Now we attempt a global search of the \(k\)-\(\beta\)-\(\rho\)-\(\eta\)-\(\mu_{IR}\) space.
We set up a design of starting points. The freeze command is described here.
freeze(
runif_design(
lower=c(Beta=5,rho=0.2,eta=0,k=1,mu_IR=0.5),
upper=c(Beta=80,rho=0.9,eta=1,k=30,mu_IR=2),
nseq=500
),
seed=2062379496
) -> guesses
We will use one of our ‘mif2d_pomp’ objects in the following. We fix \(N\) and \(\mu_IR\).
mf1 <- mifs_local[[1]]
fixed_params <- coef(measSIR,c("N"))
The following may look somewhat familar.
registerDoRNG(274481374)
foreach(guess=iter(guesses,"row"), .combine=rbind) %dopar% {
library(pomp)
library(tidyverse)
mf1 %>%
mif2(
Nmif=100,Np=2000,
params=c(unlist(guess),fixed_params)
) %>%
mif2(Nmif=100) %>%
mif2(Nmif=100) -> mf
replicate(
10,
mf %>% pfilter(Np=2000) %>% logLik()
) %>%
logmeanexp(se=TRUE) -> ll
mf %>% coef() %>% bind_rows() %>%
bind_cols(loglik=ll[1],loglik.se=ll[2])
} -> results
This global search from 500 starts took 240Â s on 250Â cpus. We append the results to our database and have a look.
read_csv("fitall_params.csv") %>%
bind_rows(results) %>%
filter(is.finite(loglik)) %>%
arrange(-loglik) %>%
write_csv("fitall_params.csv")
results %>%
filter(is.finite(loglik)) %>%
filter(loglik>max(loglik)-50) %>%
bind_rows(guesses) %>%
mutate(type=if_else(is.na(loglik),"guess","result")) %>%
arrange(type) -> all
pairs(~loglik+k+Beta+eta+rho+mu_IR, data=all, pch=16, cex=0.3,
col=ifelse(all$type=="guess",grey(0.5),"red"))
read_csv("fitall_params.csv") %>%
filter(loglik>max(loglik)-50) -> all
all %>%
filter(loglik>max(loglik)-10) %>%
ggplot(aes(x=k,y=loglik))+
geom_point()+
labs(
x=expression(k),
title="poor man's profile likelihood"
)
all %>%
filter(loglik>max(loglik)-10) %>%
ggplot(aes(x=eta,y=loglik,color=k))+
geom_point()+
labs(
x=expression(eta),
title="poor man's profile likelihood"
)
Perhaps unsurprisingly, these poor-man’s profiles suggest there is a very flat ridge in the surface. We can see some of the tradeoffs involved with a scatterplot matrix.
pairs(~loglik+k+Beta+eta+rho+mu_IR, pch=16, cex=0.3,
data=filter(all,loglik>max(loglik)-10),
col=ifelse(round(all$k,2)==10,"blue","red"))
Because this ridge appears flat in the \(\eta\)-direction, we will profile over \(\eta\) to improve our estimates. The following sets up a design of guesses, to be used as starting points for the profile computation.
read_csv("fitall_params.csv") %>%
filter(loglik>max(loglik)-10,loglik.se<1) %>%
sapply(range) -> box
freeze(
profile_design(
eta=seq(box[1,"eta"],box[2,"eta"],length=20),
lower=box[1,c("Beta","rho","k","mu_IR")],
upper=box[2,c("Beta","rho","k","mu_IR")],
nprof=15, type="runif"
),
seed=1893696051
)-> guesses
registerDoRNG(830007657)
foreach(guess=iter(guesses,"row"), .combine=rbind) %dopar% {
library(pomp)
library(tidyverse)
mf1 %>%
mif2(
Nmif=100,Np=2000,
params=c(unlist(guess),fixed_params),
rw.sd=rw.sd(Beta=0.02, rho=0.02, mu_IR=0.02, k=0.02),
partrans=parameter_trans(log=c("Beta","mu_IR","k"),logit="rho"),
paramnames=c("Beta","mu_IR","k","rho")
) %>%
mif2(Nmif=100,cooling.fraction.50=0.3) -> mf
replicate(
10,
mf %>% pfilter(Np=2000) %>% logLik()
) %>% logmeanexp(se=TRUE) -> ll
mf %>% coef() %>% bind_rows() %>%
bind_cols(loglik=ll[1],loglik.se=ll[2])
} -> results
results %>%
filter(
is.finite(loglik),
loglik.se<1
) %>%
group_by(eta) %>%
filter(rank(-loglik)<=2) %>%
ungroup() %>%
gather(variable,value,loglik,mu_IR,rho,Beta,k) %>%
ggplot(aes(x=eta,y=value))+
geom_point()+
labs(y=NULL)+
facet_wrap(~variable,scales="free_y")
Licensed under the Creative Commons Attribution-NonCommercial license. Please share and remix noncommercially, mentioning its origin.